
    ### Estimators

    # Set na.rm=TRUE to accommodate the Trim method
    Mean = function(x) mean(x,na.rm=TRUE)
    Sum  = function(x) sum(x,na.rm=TRUE)

    HT = function(N,eX,Z,Y){
        return( Mean(Z * Y / eX) - Mean((1-Z)*Y/(1-eX)) )
    }
    
    Ratio = function(N,eX,Z,Y){
        return( Sum(Z*Y/eX)/Sum(Z/eX) - Sum((1-Z)*Y/(1-eX))/Sum((1-Z)/(1-eX)) )
    }
    
    DR   = function(N,eX,Z,Y,Yhat1,Yhat0){
        part1 = Sum(( Z*Y - (Z-eX)*Yhat1 )/eX) / N
        part0 = Sum(( (1-Z)*Y + (Z-eX)*Yhat0 )/(1-eX)) / N
        return(part1-part0)
    }

    
    Imb = function(N,eX,Z,X){
        
        Z = Z[!is.na(eX)];  X = X[!is.na(eX),];    eX = eX[!is.na(eX)]
        N = length(Z)
        
        X.full = cbind(rep(1,N),X)
        
        part1 =   t(as.matrix(Z / eX - (1-Z) / (1-eX))) %*% X.full / N
        part2 =   t(X.full) %*% X.full / N
        Imb10 =   sqrt(part1 %*% solve(part2) %*% t(part1))
        
        return(Imb10)
        
    }
    
    
    
   
    
    
    
    
    
    
    